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ABSTRACT 



Two numerical models were developed which generate synthetic 
seismograms in formats similar to those of records obtained in field 
exploration. The models solve the hyperbolic wave equations in one 
and two dimensions by finite difference approximations in initially 
undeformed solution domains of transversely-isotropic layered media 
subjected to a time varying, dilatational forcing stress applied at 
the surface. A velocity attenuation term was included in the models 
to inhibit strong boundary reflections. 

The one-dimensional model produces synthetic reflection profiling 
seismograms for arbitrary horizontal or dipping layers. The two- 
dimensional model generates synthetic reflection/refraction seismograms 
for horizontal layered media with arbitrary distribution of wave 
velocities, Poisson’s ratio and density. Several sample records were 
produced for some representative velocity structures. The synthetic 
seismograms were interpreted and gross correlation was carried out 
as if they were actual field records. 
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I. INTRODUCTION 



A major problem in exploration geology is the determination of 
structural relationships for areas in which the rocks cannot be directly 
examined. This is usually a subjective process, involving the extra- 
polation of lithologic facies from a region where control is available, 
such as surface outcrops, mines or wells, across an unknown or 
inaccessible area to another point with known structure. Indications 
of subsurface structure, including magnetic, seismic and gravity 
measurements as well as interpretation of surface morphological 
features, are frequently relied upon to aid this process. 

The problem becomes particularly acute in oceanic regions where 
the only available control may be continental outcrops several hundreds 
of miles away, a few bottom grab samples and shallow sediment cores 
and a rough bottom topography. Remote measurements may then become 
the only method of initial structure determination. 

Frequent use is made" of the manner in which the ocean bottom and 
sub-bottom rocks reflect and transmit acoustic or seismic energy, not 
only for the correlation of near-shore with on shore structures, but 
also for initial determination of sub-bottom structures. When con- 
ducting seismic exploratory operations, the information obtained is 
usually presented as a two-dimensional graph (seismogram) on which 
three variables are shown. These variables are: horizontal distance, 

time and intensity of energy received. Depth information may be 
inferred from the time scale with a known or assumed vertical distri- 
bution of sonic velocity. A high level of intensity is usually 
assumed to represent a point or surface at which the physical 
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parameters of the media change more or less abruptly, reflecting energy 
or providing a surface along which energy might be propagated. 

Interpretation of the seismograms and the development of a sub- 
surface structure can be aided by the use of a mathematical model of 
the supposed structure which would produce a synthetic seismogram. 

The synthetic seismogram could be compared with the field record for 
correlation. 

This paper presents the development of two numerical models which 
produce synthetic seismograms for layered media subjected to an initial 
time varying, compressional (dilatational) stress input applied at the 
surface. The models were designed to have fairly wide variability 
while remaining simple enough to preclude the use of inordinate 
computer time and storage requirements. The one- dimensional model 
produces a reflection profiling seismogram. It is applicable to the 
case of vertically incident dilatational waves reflected and transmitted 
by arbitrary bottom profile and sub-bottom layered structure of hori- 
zontal or dipping interfaces. The two-dimensional model yields a set 
of seismograms applicable to a reflection/refraction problem where the 
receivers are distributed in a horizontal array. Both dilatational 
and shear stresses are treated. The region modeled is one of horizontal 
layered media with arbitrary distribution of shear and dilatational 
wave velocities and of density. 

No attempt was made here to include the wide variability of possible 
media which is an integral part of the large SLAM (Stress waves in 
LAyered Media) code model developed at the Illinois Institute of 
Technology (IIT) fconstantino 1969] primarily to treat the neighbor- 
hood of high energy (nuclear) explosively generated ground shocks. 
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The IIT model has the capability of simulating a wide variety of 
conditions but requires a large amount of computing time and space in 
order to produce a single trace. It is felt that the two models pre- 
sented here have a greater applicability to marine and terrestrial 
seismic exploration problems. 

These models were developed on the basis of numerical approxima- 
tions to the solutions of the wave equations in one and two dimensions 
by finite difference methods. The computer programs were written in 
the FORTRAN IV language and were run on a dual IBM 360/67 system. The 
plotting subroutines used to drive a CALCOMP 765 plotter were prepared 
by the staff of the United States Naval Postgraduate School Computer 
Facility ^ 1969J 
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II. THE SEISMIC REFLECTION PROFILING MODEL 



The problem of synthetic seismic profile generation was approached 
as one of repeated transmission of normally incident longitudinal 
(dilatational) disturbances through transversely-isotropic media. The 
seismogram was developed by changing the velocity structure in the 
solution domain between pulse transmissions, corresponding to a change 
in position of the shot point as the profiling vessel moves across an 
area of changing water depth and sub-bottom structure. The solutions 
are presented as successive traces as though from a pressure hydro- 
phone located at any depth below the shot point. 

A. THE MATHEMATICAL MODEL 

It was assumed that the wave velocity (C) is continuous throughout 
the domain and either varies linearly with depth (x) or is constant. 
Additionally, it was assumed that in transition zones, where the wave 
velocity does change, density remains constant. Under these assump- 
tions, the equation of motion in one dimension for dilatational waves 
can be written (modified from Wolf jl937] ) : 



where u is longitudinal particle displacement and A is a velocity 
attenuation parameter which is a function of x. Carrying out the 
indicated differentiation of the right side of equation (1) yields: 



E: 




( 1 ) 




( 2 ) 
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Four collateral conditions are required in order to complete a 
well-posed problem. The following conditions were selected: 
initial conditions: 
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c , 
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boundary 


conditions : 








Bl: 


(L,t> = 
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B2: 
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(6a) 






— » 




(6b) 



where L is the distance to the lower boundary. 

The initial conditions allow the solution domain to be initially 
undeformed and at rest. The first boundary condition requires that 
displacements are zero for all time at the lower boundary, which then 
becomes perfectly reflective. The second boundary condition allows a 
variable stress to be applied to the upper boundary for t ^ t^. For 
t > tj^, normal stresses are not transmitted across this boundary. It 
then becomes a free surface. 

The second term on the right side of equation (1),A^^ , represe- 
sents an attenuation factor, dependent on depth and particle velocity, 
which was used to reduce the magnitude of the disturbance in the lower 
regions of the domain so that reflections from the lower boundary would 
not mask late arriving multiple reflections from points nearer to the 
surface. The function A(x) is assumed to be zero in the region from 
which reflections are to be studied. 
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B, THE SOLUTION DOMAIN 



The solution domain (Fig, 1) is divided by N equally spaced mesh 
points into N-1 subintervals each of length Sx. A sample wave 
velocity distribution is shown in the figure. The upper boundary (x = 

0) is shown on the left. In the region O^x^x^, the wave velocity is 
constant and equal to Cl. For x^^. x ^ L, the velocity is equal to C2. 

In the transition zone, ^ ^ ^2* velocity varies linearly from 

Cl to C2. 

The time-depth mesh is an array of N by M mesh points. Each 
rectangular mesh area is then of size S x by This mesh is the 

full solution domain and the objective of the model is to find a solu- 
tion of equation (2) at each mesh point under the collateral conditions, 
equations (3) through (6). 

C. THE DIFFERENCE EQUATIONS 

A central difference form for each of the partial derivatives in 
equation (2) can be written in terms of the (i,j) mesh point [Smith 
1965] . Figure 2 shows the (i,j) and neighboring mesh points. The 

required derivatives, in central difference form, become: 



Ft “ 
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Figure 1. The solution domain 
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J+1 




Figure 2. Mesh detail. 
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Substituting these into equation ( 2 ) gives the difference equation of 



motion for the (i,j) mesh point: 



) ^-1 - j ^ Ot , )- 1 — C" ( Ol-\, j 






(7) 



The collateral conditions become: 
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= F(t) , 0 4 t ^ t^ , 








Fj - 0 , t > . 






Solving 


equation (7) for u^ gives: 






E: 










+ j;-i,j (c:o-3 aCcc^i-a-OR^) 
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where : 


- = (f.V • 






Writing 


equation (12) for j=l and substituting 


u^ Q from equation 


(9) 


yields 


the initial conditions: 
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12: 


- o • 




(13) 
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Writing equation (12) for i=l and substituting u , from equation (11) 

® > J 



yields the boundary conditions: 



O. 



- o 



Bl: 

(i*4A, -t) + Fjb<n^)]/[^-: 



(10a) 



(14) 



The equations (12), (8), (13), (10a) and (14) then form a well- 
posed problem for numerical solution. Initially, the j*=l and j=2 rows 
are set to zero following equations (8) and (13). The j+1 row can then 
be solved utilizing equation (10a) for the lower boundary (i=N) , 

(14) for the upper boundary (i=l) and (12) for the interior points. 

Note that the value of u anywhere along the interior of the j+1 row is 
dependent only on the values of u at three points along the j row and 
one on the j-1 row (shown as circles in Fig. 2). 



D. STABILITY AND CONVERGENCE 
1. Stability 

Because of the large number of arithmetic operations involved 
9 

(approximately 10 for the solution spaces used), it is necessary to 
insure that there is a stable decay of errors introduced in the various 
computations. Studies of the stability of the centered difference 
formulation of the basic wave equation: 






(15) 



liavo been carried out by several investigators [^Smith 1965] . Smith 

investigated the stability by the Fourier series method and showed 
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that the centered difference formulation is stable for the mesh 
ratio : 



Sv 



4 



c 



( 16 ) 



The addition of the velocity gradient and velocity attenuation 
terms in equation (2) alters the stability conditions somewhat, A 
large velocity gradient, approaching discontinuous velocity change, 
generates instability in the gradient region. Investigation of the 
stability conditions showed that if the condition: 



^c\. - 4 Ci. , (17) 

were imposed on the gradient regions, the requirement on the mesh ratio 
[equation (16)] becomes: 



St / 0-89 

^ C 



(18) 



The addition of a negative attenuation term tends to decrease 
the right side of equation (18). This would lead to a more stringent 
requirement on the mesh ratio. This effect was found to be quite small 
for the values of A(x) necessary to damp lower boundary reflections. 
Empirically, it was found that stability would be assured by requiring 
that: 



St / 0-35 

C 



(19) 



2 . Convergence 

A finite difference scheme is said to be convergent if the 
exact solution of the difference equations approaches the exact 
solution of the partial differential equations as Sx and St 
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approach zero together. Forsythe and Wasow [l967] have shown that 
convergence of a finite difference scheme such as equation (12) with 
boundary conditions such as equations (8), (13), (10a) and (14) is 
convergent as long as the mesh ratio is such that the (i, j+1) mesh 
point lies within the intersection of the characteristics through the 
(i-l,j) and(i+l,j) mesh points (see Fig, 2). The equation defining 
the characteristics for equation (2) is; 



This leads to: 



^ - -i- r 

dt 



( 21 ) 



Convergence is then assured for: 



^ / JL 

^ C 



(16) 



It is apparent then that as long as equation (19) holds, the finite 
difference formulation will be both stable and convergent. 



E. DESCRIPTION OF THE PROGRAM 

The program requires 140K bytes of main core storage. An average 
run with a 251x1800 mesh grid with 41 shot points requires approxi- 
mately 13 min of main run time exclusive of output time required for 
the plotting subroutines. 

1 . Distribution parameters 

The advantage of this model lies in the wide variability of the 
conditions which may be easily simulated. The primary parameters used 
to model various conditions are the distribution parameters. They are: 
a. wave velocity distribution, C(x,t), 
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h. forcing function, F(t), 

c. attenuation parameter, A(x). 

The forcing function and the attenuation parameter are put intt> 
tile model as s ingle-valued , continuous, polygonal line functions of 
tile indicated dependent variable. The wave velocity is put into the 
model as a two-dimensional array of points which describe the depths to 
isovelocity layers at successive horizontal stations. Figure 3 shows a 
sample distribution of isovelocity horizons with horizontal stations 
and shot points indicated. Horizontal station distances and the 
velocity defining array for this sample are shown in Tables I and II, 

At each shot point, the wave velocity is a single-valued, continuous 
polygonal line function of depth. An example of wave velocity distri- 
bution at a single shot point is shown in Fig. 4 along with a sample 
attenuation parameter and forcing function. The polygons for A(x) and 
F(t) are defined in Tables III and IV. 

2 . Numerical parameters 

The numerical parameters, which remain constant throughout the 
problem, are used to control the computation of the distribution para- 
meters and to set the arrangement and size of the solution domain. 

They are: 

a. mesh size in the x direction, (DX) , 

b. mesh size in the t direction, (DT) , 

c. horizontal distance between successive solution stations 
(shot point spacing), (DK) , 

d. number of total x mesh points, (N), 

e. number of total t mesh points, (M) , 

f. number of solution stations, (KK), 
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Figure 3. Sample velocity structure for seismic profiling 
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Figure 4, Sample distribution parameters 

(a) Wave velocity at one shot point. 
(h) Attenuation parameter, 

(c) Forcing function. 
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g. number of points defining the distribution parameter 
polygons : 

(1) wave velocity polygon, (NC) , 

(2) attenuation parameter polygon, (NA) , 

(3) forcing function polygon, (NF) , 

h. number of horizontal stations in the velocity defining 
array, (NG) , 

i. mesh point at which the output stress is to be measured, 

(IP), 

j. horizontal distance to the first shot point, (EEE) . 
Additional numerical parameters are put into the program which 

control the size, scale and arrangement of the output traces. For 
details of their use see the explanatory notes in the program. 



Table I. 



Sample horizontal 

Station 

1 

2 

3 

4 



station distances 

Distance 

0.0 

2.0 

10.5 

20.0 



Table II. 

Sample velocity defining array 



Station Numbers 





1 


2 


3 


4 


Velocity 


Depth 


0.0 


0.0 


0.0 


0.0 


Cl 


to 


2.0 


2.0 


4.5 


9.8 


Cl 


isovelocity 


2.2 


2.2 


4.7 


10.0 


C2 


horizons 


4.7 


4.7 


4.7 


10.0 


C2 




4.9 


4.9 


4.9 


10.2 


C3 




20.0 


20.0 


20.0 


20.0 


C3 
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Table III. 



Sample attenuation parameter polygon 



Depth 


A(x) 


0.0 


0.0 


10.5 


0.0 


20.0 


-2.0 



Table IV 



Sample forcing function polygon 



Time 


F(t) 


0.0 


0.0 


0.10 


0.0 


0.40 


-0.15 


0.70 


0.0 


2.0 


0.0 



F. SCALING THE PROBLEM 

All wave velocities and depths were put into the model in non- 
dimensional form. The wave velocity in the upper layer and the depth 
to the first transition zone were chosen as the definitive values for 
non-dimensionalizing the problem. Non-dimensional time was then 
determined by these values through the mesh ratio. Practically, the 
scale of the problem is determined by the wave velocity values and one 
of the mesh sizes. In general, the wave velocity is entered as C 
distance units per time unit. For a mesh ratio of: 



St _ O - 025 _ i 

O-iO ^ 



( 22 ) 



the stability condition [equation (19)] allows the wave velocity to 
vary from 0.0 to 3.4 distance units per time unit. If a value of C of 
0.8 is to represent dilatational wave velocity in water (1.5 km/sec), 
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then the scale of the problem would depend on a choice of either time 



or depth mesh size, chosen in view of the definition desired, 
shows some scales where the time mesh size was specified. 



Table V 
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Numerical parameters: N = 251, M = 2000 
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III. THE REFLECTIQN/REFRACTION MODEL 



Synthetic reflection/refraction seismogram generation requires 
solution of the wave equation in at least two dimensions. A rectangu- 
lar solution space was chosen with three rigid, reflective boundaries 
and one free boundary corresponding to the air-water or air-ground 
interface. In order to produce a model with wide variability as well 
as one that may be easily manipulated without inordinate computer time 
and storage requirements, it was assumed that the media in the solution 
space is horizontally isotropic. The y-derivatives of density and 
Lame*s constants (and hence of shear and dilatational wave velocities) 



are continuous throughout the solution domain. Further, it was 
assumed that density and the dilatational and shear wave velocities, 

C and S, are either constant or vary linearly with depth. Under these 
assumptions, the equations of motion in two dimensions for the propaga- 
tion of shear and dilatational disturbances can be written (modified 



will then be zero 



A. THE MATHEMATICAL MODEL 



It was assumed that density, p , and Lame*s constants, and 



from Officer 




and Takeuchi 




): 



El: 





(23) 



E2: 




( 24 ) 
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where u and v are the particle displacements in the x and y directions 



and B is a velocity attenuation factor. The stresses are: 



Dilatational stress : 



Vc)v ; 



Shear stress: 



/ V5>! 



These stresses are computed in the program and are used, in addition 
to displacements, to form the seismograms. 

Carrying out the indicated differentiation on the right side of 
equations (23) and (24) and collecting terms yields: 
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ct^ 


= ( 


b +v<^ 

P' ^ 




, i 


bb 


/bt) , 










bH ) 




_ / 


^ 'l/i 






' \ 




y 


+ 1- 








P 






bM j 



"u 









( 






(25) 
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where it has been assumed that: 






^ M ^ o ? 



- o/t _ 



^ "o 



_ 



o 









_ 

b-i bx 



31 



The eight collateral conditions chosen to complete a well-posed 



problem were: 



initial conditions : 
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(X, 1-1 .o') 


= O 


(27) 


12: 
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(28) 


13: 
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boundary 


conditions : 
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(33) 




where: F(y,t) = 


0 for t t^ 




B4: 


+ Vy 


^ o 


(34) 



The initial conditions allow the solution space to be initially unde- 
formed and at rest. The first two boundary conditions (Bl and B2) 
require that displacements are zero for all time at the subsurface 
boundaries, which then become perfectly reflective. The other boundary 
conditions (B3 and B4) allow the x=0 boundary to become a free surface. 
Shear and dilatational stresses are not transmitted across the surface 
for all time except for an externally applied, time-varying dilatational 
stress, applied at a point or combination of points for t ^ t^. 

As in the one-dimensional model, an attenuation term was included 
in the equations of motion (El and E2) in order to keep boundary 
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reflections from masking multiple reflections in the central regions. 
The function B(x,y) is non-zero only in the vicinity of the three rigid 
boundaries and is zero in areas where disturbances are to be studied. 

B. THE SOLUTION DOMAIN 

The solution domain is a three-dimensional rectangular array of L 
by M by N mesh points dividing the domain into (L-1) by |M-1) by (N-1) 
equal subdivisions, each of size Sx by Sy by ^t. Figure 5 shows the 
x-y solution space. The total solution domain is formed from N succes- 
sive solution spaces or planes. A detail of the distance-depth- time 
mesh in the vicinity of mesh point (ijj^k) is shown in Fig. 6. 

C. THE DIFFERENCE EQUATIONS 

As in the one-dimensional model, a centered difference form for 
each partial derivative in equations (25^ and (26) can be written. 
Making use of the definitions for wave velocities: 



finite difference forms of equations (25) and (26) may be written in 
terms of wave velocities and their derivatives. Carrying out this 
substitution and solving for the time leading (i,j,k+l) term for u and 
V yields: 




(35) 
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( 37 ) 
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Figure 6. Isometric detail of three- 
dimensional mesh. 
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( 38 ) 
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A well-posed problem for numerical solution is then formed by 
equations (37) and (38) with collateral conditions (39) through (46). 
Initially, the k=l and k=2 planes are set to zero following equations 
(39) through (42), The k+1 plane can then be solved using equations 
(43) through (46) for the boundaries and equations (37) and (38) for 
the interior mesh points. 



D. STABILITY AND CONVERGENCE 
1. Stability 

Under the conditions that the attenuation and velocity gradient 
terms are zero and that Sx = Sy, the difference equations (37) and 
(38) reduce to those for which a stability analysis was carried out by 
Alterman and Karal [l968] , and Alterman and Rotenberg [l969] 

The former paper showed that the stability criterion is given by: 

^ ^ i. 

8X c 



1 ^ 



(ir 



(47) 



In the latter paper, th*e authors indicated that stability was guaranteed 
under the condition that: 



/ O * Q(d 

with the assumption that S/C = 0.55, 
of 0.28. Poisson’s ratio is defined 



cr - 



i (1)^- 

cr-i 



’ (48) 

corresponding to a Poisson’s ratio 
in terms of wave velocities as: 

(49) 



Empirically, it was found that the addition of attenuation and 
velocity gradient terms reduced the right side of equation (47), while 
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an increase in 



Sy ( Sy > Sx) reduced the left side. An approximate 



stability criterion was established, given by: 




(50) 



under the condition that the restrictions on velocity gradients (for C 
and S) given by equation (17) apply. 

2 . Convergence 



long as the mesh ratios were sufficiently small to insure that the 
(i,j,k+l) mesh point lies within the characteristic surfaces through 
adjacent mesh points in the k plane. Fox equated the characteristic 
surfaces with successive positions of a wave front. The most restrictive 
condition on convergence is then given by the maximum dilatational 
wave velocity in the solution space. This would give a less stringent 
mesh ratio condition than the stability condition. Therefore, as 
long as equation (49) holds, the method will be convergent as well as 
stable. 

E. DESCRIPTION OF THE PROGRAM 

The program produces a set of four seismograms: shear and dilata- 
tional stress and vertical and horizontal displacement. The data is 
written on a temporary storage device (tape, disk or drum). A second 
program then reads the data and plots the seismograms. An average run 
with a 71x81 grid, 800 time mesh points and 25 detectors requires 278K 
bytes of main core storage and a run time of 18 min. The output 
requires 325K bytes of temporary storage space. The primary limiting 



Following the methods of Forsythe and Wasow 




and of 



Fox et 




, it was assumed that convergence would be assured as 
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factor on a given problem is the time required for input/output to and 
from the temporary storage device. 

1. Distribution Parameters 

As in the one- dimensional model, the distribution parameters 
are used to model the physical conditions in the solution space. They 
are: 

a. dilatational wave velocity, C(x), 

b. shear wave velocity, S(x), 

c. density, RH(x) , 

d. forcing function, F(y,t), 

e. attenuation parameter, B(x,y). 

Wave velocities and density are put in the program as single- 
valued, continuous, polygonal line functions of depth. The forcing 
function is put in as was done in the one-dimensional model along with 
horizontal stations at which F is to be applied. The attenuation 
parameter is put into the program as two single-valued, continuous, 
polygonal parameters. They are: 

el. vertical attenuation parameter, BA(x), 

e2. horizontal attenuation parameter, BB(y). 

The full attenuation parameter, B(x,y) , is then computed as follows: 

= MAX.[BA(x), - (51) 

Tables VI through VIII give sample numerical input used to 
describe the various distribution parameters. Their distributions in 
the solution domain are shown graphically in Fig. 7 through 9. 
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Table VI 



Sample wave velocity and density polygons 



Depth 


C(x) 


S(x) 


RH(x) 


0.0 


0.6 


0.0 


1.0 


3.0 


0.6 


0.0 


1.0 


3.3 


1.2 


0.82 


2.1 


7.0 


1.2 


0.82 


2.1 



Table VII. 



Sample forcing function polygon 



Time 


F(t) 


0.0 


0.0 


0.1 


0.0 


0.7 


-0.15 


1.30 


0.0 


2.50 


0.0 



Table VIII. 



Sample attenuation parameter polygons 



Depth 


BA(x) 


Distance 


BB(y) 


0.0 


0.0 


0.0 


-2.5 


4.0 


0.0 


1.0 


-1.0 


7.0 


-2.5 


2.0 


0.0 






10.0 


0.0 






11.0 


-1.0 






12.0 


-2.5 



40 




41 



Figure 7, Sample distribution of wave velocities and density. 
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Sample distribution of the attenuation paramc^ter. 




2 , Numerical Parameters 



These parameters are used to control the computation of the 
distribution parameters and to set the arrangement and size of the 
solution domain. They are: 

a. mesh size in the x direction, (DX) , 

b. mesh size in the y direction, (DY) , 

c. mesh size in the t direction, (DT) , 

d. number of total x mesh points, (L) , 

e. number of total y mesh points, (M) , 

f. number of total t mesh points, (N), 

g. number of points defining the distribution parameter 
polygons ; 

(1) wave velocity and density polygons, (NC), 

(2) forcing function polygon, (NF) , 

(3) attenuation parameter polygons: 

(a) vertical, (NBA), 

(b) horizontal, (NBB) , 

h. X mesh point at which the output is measured: 

(1) stress, (lA), 

(2) displacement, (ID), 

i. number of horizontal output stations, (KK) , 

j. y index of the first output station, (JB) , 

k. mesh increment between output stations, (JJ), 

l. shot point or area location: 

(1) first shot mesh point, (MA) , 

(2) last shot mesh point, (MB), 

m. last time mesh point at which F(y,t) is to be applied, (LL) . 
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F. SCALING THE PROBLEM 



The scale of the problem is determined by choices of wave velocity 
values and either &t or 8x and Sy- For mesh ratios of: 



St _ Q OC ^ ^ 

Sx. O.iO 2. 



(52a) 



St - Q-Q9 - J_ 

O' ^ 



(52b) 



and under the condition on shear velocity that S = C/ V3" (CT - 0.25), 
the dilatational wave velocity may vary from 0.0 to 1.86 distance units 
per time unit, and shear wave velocity may not exceed 1.07 distance 
units per time unit. If a value of C = 0.60 is to represent dilata- 
tional wave velocity in water of 1.5 km/sec, then the scale of the 
problem would depend upon a choice of either time or distance mesh 
size. Table IX shows scales available for various distance mesh sizes 
specified. 
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Table IX 
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IV. INTERPRETATION 



The synthetic seismic records produced by the two models presented 
in this paper were interpreted first in the manner of actual field 
records without direct reference to the structural profiles used. 
Extensive use was made of gross seismic interpretation techniques pre- 
sented by Grant and West [l965] , Dobrin [l960] , and Bullen [l959] . 

Interpretation of the reflection and transmission of incident pulses 
was aided by the results of Alterman and Karal [l96S] , Alterman and 
Rotenberg [l969j , Gupta [ 1966a and 1966b] , Sengbush, Lawrence and 
McDonal [l96l] and Tooley, Spencer and Sagoci [l965j . 

A. THE REFLECTION PROFILING SEISMOGRAMS 

Each reflection profile produced by the model is designated by an 
alpha-numeric sequence. For example, C3H F7 A14 identifies the exact 
structure (C3H), forcing function (F7) and attenuation parameter (A14) 
used. 

Six profiles were developed using three basic structures. The 
structures, shown in Fig. 10 through 15, can be described in terms of 
marine geology as follows: 

(1) Shelf break with horizontal layers. 

(a) Increasing velocity in each successive layer (C3H). 

(b) With one low velocity layer (C4D). 

(2) Flat bottom with double wedge. 

(a) High velocity wedge (C5A) . 

(b) Low velocity wedge (C6C). 

(3) Flat bottom, horizontal layering, penetrated unconformably by 
a dipping bed. 
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(a) High velocity bed (C7A) . 

(b) Low velocity bed (C8A) . 

All traces utilized forcing function F7 (Fig, 16) and attenua- 
tion parameter A14 (Fig. 17). 

1 . Variable Depth Records 

As an aid to interpretation, several single traces were 
developed with a single velocity structure, where the depth of the 
detection point was varied between traces. The traces are shown in 
Fig. 18 along with the mesh point at which each was computed. These 
various detection points are shown as circles in Fig. 10. The traces 
were generated using a source located at a horizontal distance of 0.5 
units for the problem designated C3H F7 A14. The first disturbance in 
each trace, designated d, is the direct wave. and designate the 

reflected pulses from the first and second transition zones, while the 
first primary multiple (P^ reflected from the free surface) is shown 
by Other multiples are marked in a similar manner. 

These traces show’clearly the development of a train of alter- 
nating negative and positive pulses behind the main pulse. The change 
in pulse shape is also quite evident. As the pulse progresses, the 
amplitude decreases and the pulse width increases. Additionally, the 
base and tip of the pulse form tends to become more rounded. These are 
primarily mesh effects and are dependent on mesh size (DX) . When one 
particle is unloaded as the pulse passes by, the adjacent particle is 
still under stress, causing a reverse stress between the particles. In 
gradient zones, where the derivative of C is non- zero, there is also an 
elastic rebound effect which reinforces the pulse train. 
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Figure 16. Forcing function F7 
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Figure 17. Attenuation parameter A14 
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2 . The reflection profiling seismograms 



Synthetic reflection profiling seismograms are shown in Fig. 19 
through 26. Figure 19 is an uncorrelated version of Fig. 20, included 
so that lineations may be observed without the influence of correlation 
lines shown on the remaining figures. Solid lines on the seismograms 
mark lineations which were correlated with first reflections from a 
particular horizon. Dashed lines indicate multiple reflections which, 
if used as primary reflections, might lead to erroneous interpretations 
of structure. For each correlation, a particular phase was chosen and 
followed laterally through successive traces. Figures 20 and 21 show 
clearly the false dips introduced by the sloping bottom, increasing the 
pulse travel time in the lower velocity upper layer. The bottom slope 
causes a change in thickness of the top two layers (water and the first 
subbottom layer) . The false dip introduced is a function of the bottom 
slope and the ratio of velocities involved. Note that there is a change 
in dip at each point where a layer pinches out. 

Figures 22 and 23 show similar false dips and horizons intro- 
duced by variations in the thickness of a single layer. The high 
velocity wedge causes a doming up to be introduced into a horizontal 
transition zone below, while the low velocity wedge has the reverse 
effect. The most obvious feature of the records for velocity structures 
C7A and C8A (Fig. 24 and 25) is the relatively high amplitudes in the 
records in the areas where the dipping beds pinch out against the flat 
bottom. This is caused by a combination of factors. In velocity 
structure C7A (Fig. 24), there is a sharp change in velocity contrast 
across the bottom as the dipping bed outcrops. The steeper gradient 
then reflects more energy, causing an increase in amplitude in both 
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Figure 19 

Uncorrelated reflection profile C3H F7 A14 
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Figure 19. Uncorrelated reflection profile C3H F7 A14 



Figure 20 

Reflection profile C3H F7 A14 



DX: 


0.10 


DT; 


0.025 


DK: 


0.50 


N: 


251 


M: 


1800 


IP: 


2 



Horizontal distance (d) 
versus 
time (t) 

First reflections 
Multiple reflections 



57 




58 



Figure 21 

Reflection profile C4D F7 A14 
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Figure 22 



Reflection profile C5B F7 A14 
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Figure 23 



Reflection profile C6C F7 A14 
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Figure 24 

Reflection profile C7A F7 AlA 
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Figure 25 

Reflection profile C8A F7 A14 
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Figure 26 

Precision reflection profile C4D F7 Al4 
DX: 0.10 DT: 0.025 DK: 0.10 FEE: 4.0 

N: 251 M: 1800 IP: 2 
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first reflections and multiples. Additionally, there is a resonance 
effect of the large number of multiples being generated. This resonance 
effect is also present in velocity structure C8A (Fig. 25). There is 
also a tuning effect due to the relationship between the thickness of 
the transition zone and the primary period of the incident pulse 
Sengbush, Lawrence and McDonal 196l] which would lead to high ampli- 
tude reflected pulses. The possibility that instabilities were intro- 
duced in these high gradient areas was considered. All parameters, 
however, are well within stability criteria established in Section II. 

Figure 26 is an example of a "precision trace^’ that may be 
produced by the model. This seismogram presents a section of a record 
obtained with the same velocity structure as that shown in Fig. 21 
(C4D F7 A14) with an expanded horizontal scale. The horizontal 
distance covered includes the region where the first sub-bottom layer 
pinches out. This area (from 4.0 to 8.0) is marked on Fig. 21 with a 
horizontal line in the record. The "precision trace" was obtained by 
reducing the shot point spacing (DK) . 

B. THE REFLECTION/REFRACTION SEISMOGRAMS 

Each set of synthetic seismograms developed by the model is 
designated by a letter-number sequence similar to that used in the 
reflection profiling model. CA F5A B3 L2 identifies the exact wave 
velocity and density structure (CA) , forcing function (F5A) , attenua- 
tion parameter (B3) and shot point location (L2) used to generate the 
records. The seismograms presented in this section were developed 
using the following velocity structures (Fig. 27 through 29): 

(1) Two elastic layers, (CA) . 
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o o Refraction spread 



Figure 27. Velocity structure CA 




i 9 Bottom spread 

Figure 28. Velocity structure CB 
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Figure 29. Velocity structure CC 
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(2) Water layer over elastic layer, (CB) . 

(3) Three elastic layers, (CC). 

With the exception of the water layer in CB, all velocity structures 
were designed with Poisson* s ratio of 0,25. 

All traces were produced utilizing forcing function F5A (Fig. 30), 
attenuation parameter B3 (Fig. 31> and shot point location L2. At this 
shot point location, the surface dilatational forcing stress is applied 
at mesh points j=15 and j=16. This is marked by the letter *’F” in 
Fig. 27 through 29. 

1. The Seismograms 

Selected correlated synthetic seismograms are shown in Fig. 32 
through 36. Figure 32 is a complete set of four S 3 mthetic seismic 
records for a close (reflection) spread of 25 surface detectors for 
velocity structure CA. The detector spread is shown in Fig. 27. Com- 
parison of these records shows the interdependence of the displacements 
and stresses. Actual land field records would include only vertical 
and/or horizontal displacements, or particle velocity, while marine 
records would only consist of the dilatational (pressure) record. The 
additional synthetic seismograms enable the proper phase arrivals to be 
quickly identified. A set of theoretical time-distance curves was 
developed for velocity structures CA and CB (Fig. 37), where it was 
assumed that the width of the transition zone was zero. There is good 
agreement between these theoretical arrival times and correlated linea- 
tions on the synthetic records. 

The first trace on the dilatation record corresponds to a 
normally incident wave. It bears a very close resemblance to traces 
obtained with the one-dimensional model for similar velocity structures. 
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Figure 30. Forcing function F5A 
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Figure 31. Attenuation parameter B3 
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Figure 32 



Reflection seismogram CA F5A B3 L2 
A; Dilatational stress 
B; Shear stress 
C: Vertical displacement 

D: Horizontal displacement 
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Figure 32B, Reflection Seismogram CA F5A B3 L2 : Shear stress 
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Figure 32D; Reflection Seismogram CA F5A B3 L2; Horizontal Displacement 
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Figure 33 

Refraction seismogram CA F5A B3 L2 
Shear stress 

0.10 DY: 0.15 DT: 0.05 

71 M: 81 N: 800 

2 ID: 1 JB: 22 JJ: 2 

Horizontal distance (d) 
versus 
time (t) 
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First reflections 
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Figure 34 

Reflection seismogram CB F5A B3 L2 
Dilatational stress 

DY: 0.15 DT: 0.05 

M: 81 N: 800 

ID: 1 JB: 16 JJ: 1 

Horizontal distance (d) 
versus 
time (t) 

Direct wave arrivals 
First reflections 
Multiple reflections 
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Figure 35 

Bottom detector seismogram CB F5A B3 L2 
Dilatational stress 

DX: 0.10 DY: 0.15 DT: 0.05 

L: 71 M: 81 N: 800 

lA: 20 ID: 20 JB: 16 JJ: 2 

Horizontal distance (d) 
versus 
time (t) 

0 Direct wave arrivals 

Surface reflected multiples 

Interface wave conversion 
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Figure 36 

Reflection seismogram CC F5A B3 L2 
Dilatational stress 
DX: 0.10 DY: 0.15 DT: 0.05 

L: 71 M: 81 N: 800 

lA: 2 ID: 1 JB: 16 JJ: 1 

Horizontal distance (d) 
versus 
time (t) 

0 Direct wave arrivals 

First reflections 

Multiple reflections 

Shear wave conversion 
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® Dilatational arrivals (CA and CB) 




Figure 37, Theoretical time-distance curves for 
velocity structures CA and CB, 
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Elastic rebound is present, causing the development of a cyclic wave 
train. Phase reversal is observed for events reflected from the free 
surface on all records. These features are particularly apparent in 
the contour plots discussed in the next section. 

Although no refracted arrivals are shown on any of the records 
of Fig. 32, a close examination of the traces revealed their presence 
after a refraction record was obtained. Figure 33 is the shear stress 
refraction record for velocity structure CA, The surface detector 
spread (shown in Fig. 27) extended across most of the solution space 
into the attenuation region on the right. 

In addition to the arrival of the refracted shear wave, Fig. 33 
shows the development of a small shear stress in advance of the first 
shear arrival. This is generated by the passage of the dilatational 
wave across the mesh, a result of the coarseness of the mesh relative 
to the wave velocity. As the dilatational wave passes a mesh point at 
an oblique angle relative to the x or y axis, shear strains are set up 
between mesh points. In fact, it is just this effect which causes the 
generation of the initial shear wave at the forcing function applica- 
tion point, or, more precisely, between one application point (where 
the applied dilatational stress causes some displacement') and the 
neighboring mesh point where displacements are different. The shear 
wave is then generated from two centers. For shot point location L2, 
these centers are (fictitious) mesh points (%, 14^) and(%, 16%). 

Figure 34 is the dilatational stress reflection record for 
velocity structure CB. This structure (shown in Fig. 28) is identical 
to CA except that, in the upper layer, shear wave velocity (and hence. 
Lame's constant >6^) was set to zero. The record is, of course, very 



90 



similar to that obtained from structure CA except that the amplitudes 
of the first arrivals are somewhat greater. This is caused by the lack 
of rigidity in the upper layer leading to a loss of the attenuation 
effect of shear strains on dilatational strains. 

An example of a seismogram obtained with a spread of detectors 
on the bottom of a water layer is shown in Fig. 35. Because of the 
relationship between shot point location and detector location, the 
only lineations available for correlation, except for the direct wave, 
are multiple reflections between the bottom and the surface, and re- 
fracted waves. An additional effect observed on this record is the 
generation of a head wave in advance of the direct wave. This wave is 
generated at the interface between the two layers by the dilatational 
wave travelling through the lower medium at a high velocity. 

A three-layer reflection seismogram is shown in Fig. 36. The 
main feature of this record, not seen in the two-layer seismogram, is 
the mutual interference between the first multiple from the first 
transition zone and the reflected wave from the second transition zone. 
As these two waves are out of phase with one another, the initial 
effect is one of destructive interference. As distance from the shot 
point is increased, the waves arrive from somewhat different directions, 
altering the phase relationships, leading eventually to constructive 
interference . 

2 . The Contour Output 

An on-line printer contour plotting subroutine was developed so 
that the reflection and transmission of elastic waves might be observed 
in a qualitative manner. At chosen intervals, the values of horizontal 
or vertical displacement, dilatational or shear stress or vorticity are 
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computed for each mesh point, converted to an alphabetic representation 
and printed as a rectangular contour plot. Figures 38 and 39 show 
portions of such plots, dilatational stress plotted for the problem 
CB F5A B4 L2 at time mesh levels K=16 and 91 respectively. Figure 38 
shows the generation of the initial dilatational wave by the forcing 
function. Figure 39 shows the progress of the dilatational wave as it 
is reflected and transmitted at the transition zone. Some minor modifi- 
cation of the main two-dimensional program was necessary to produce this 
form of output. Examination of the subroutine shows what additional 
parameters were necessary. Table X shows the scaled values represented 
by the plotted characters. 

An animated motion picture was made of dilatational and shear 
stress for this problem. This presentation enables the progress of 
wave fronts throughout the solution space to be readily observed, 
including phase changes where they occur. This output allowed the 
efficacy of the attenuation parameter to be determined quite easily 
during the development phase of the model. A copy of this motion 
picture on Super 8mm film is on file at the Naval Postgraduate School, 
Department of Oceanography (Code 58Ad), and is available to interested 
investigators on a loan basis. 

Some effects were noted on the contour outputs before they were 
identified on the seismograms. The dilatational plots first showed the 
development of a head wave in the upper layer (Fig. 40^ . Other pre- 
sentations clearly show the development of a shear wave in the lower 

medium due to the incidence of a dilatational wave against an elastic 
transition zone (Fig. 41) and the generation of an interface wave 

travelling along the fluid-elastic boundary in advance of the shear wave 
(Fig. 42). 



92 



DILATATION: CB F5A B4 L2 
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Figure 38. Portion of the printer contour plot showing 
generation of the dilatational wave. 

(see Table ,X for scaled values) 
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Figure 39. Portion of the printer contour plot 

showing transmission of the dilatational 
wave into the lower media 

(see Table X for scaled values) 
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Table X. 
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-0.0068750 
-0.0064583 
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-0.0056250 
-0. 0052083 
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Figure 40. Portion of the printer contour plot 
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Figure 41, Portion of the printer contour plot 

showing generation of a shear wave at 
the transition zone by the dilatational 
wave. 

(see Table X for scaled values) 
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Figure 42. Portion of the printer contour plot 
showing development of an interface 
shear wave. 

(see Table X for scaled values) 
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V . SUMMARY 



The two seismic models presented here were developed primarily as 
tools to assist in marine geophysical field exploration. The models 
produce synthetic seismograms similar in format to that of actual field 
records in order to facilitate direct comparison. Although both models 
present two-dimensional records (horizontal distance versus time) , the 
reflection profiling model is made up of a series of one-dimensional 
model solutions. The complete solution simulates a moving source of 
vertically incident pulses where the velocity structure is slightly 
altered between transmissions. In contrast, the reflection/refraction 
model treats the propagation of energy in two dimensions. Each complete 
solution simulates a single pulse travelling through the solution space. 
The synthetic seismograms are produced as records from a series of 
detectors in a horizontal array. 

The models are based on, numerical solutions to the wave equations 
in one and two dimensions in bounded domains. The equations are 
expressed in unsimplified forms which are not amenable to direct 
analytical solution. The stability of these numerical solutions and 
their convergence to the analytical solutions was investigated. Sta- 
bility and convergence criteria were established in terms of mesh 
ratios and wave velocities. 

In order to provide a wide variability of physical simulation, the 
solution domains are modeled by a number of easily varied non-dimen- 
sional parameters. These parameters are used to set the size and sub- 
division of the solution space as well as the distribution of elastic 
parameters and the shape and extent of the forcing function. 
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The structures modeled in this paper were, in general, based upon 
realistic distributions of the non-dimensional elastic parameters in 
marine environments. Some of these models have inherent resonance with 
the type and width of the forcing function used and the mesh sizes. 
Records are produced with some features which reflect this resonance. 

In particular, there is some tuning of the pulses in the narrow transi- 
tion zone in some of the models. Care should be taken when trying to 
extrapolate features found in these records to cover a more general 
case. 
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VI. FUTURE WORK 



The applicability of these models to a wide variety of studies was 
the primary consideration in their design. The interpretation of various 
velocity structures in this study was limited to a gross overview of 
simpler problems that might be encountered in field exploration. Al- 
though a full presentation of the use of all of the variable parameters 
was not attempted here, there are areas in both models where improvements 
could be made which would increase their variability. 

A first step in improvement of the one-dimensional model might be 
the inclusion of a variable depth source. It was felt that this was not 
warranted in the model due to the complexity of a forcing function at 
depth that could not be included as part of one boundary condition. 

The same improvement might be applied to the two-dimensional model. 
However, this would involve an additional order of magnitude of com- 
plexity and the model could no longer be called ’’simplified.” A more 
reasonable improvement could be made which would remove the requirement 
that all transition zones be horizontal. The wave velocities and 
density could be put into the model as two-dimensional arrays instead 
of vectors. Terms involving the y-derivative of Lame's constants and 
of density would then be retained. This improvement would approximately 
double the main core storage required (to 550K bytes) and increase the 
main program run time by about one-third (to 24 min). 

Although the attenuation parameters used in both models are con- 
sidered adequate for the purpose for which they were included (elimi- 
nation of boundary reflections), they are by no means definitive. 
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These models could be used to study the effects of various velocity 
attenuation functions by entering A(x) or B(x,y) as functions of other 
physical parameters (e*g. overburden pressure, Lame’s constants or 
particle displacement ) • 

Only a single forcing function is presented for each model. Several 
functions were examined in the development of the models including a 
flat-topped pulse, spike pulse and ’’sine” function. The models could 
be used to study the reflection and/or refraction of a wide variety of 
forcing functions, limited only by the size of the time mesh increment 
(DT). 

In Section III, it was noted that an interface wave was observed in 
shear stress on the contour plot output. This model could certainly be 
used to study the various classes of such waves under a variety of 
conditions. Some preliminary work was done in generating particle 
motion plots from the vertical and horizontal displacement outputs 
available. Although elliptical particle motions were observed as pre- 
dicted [Alterman and Rotenberg 1969] , limitations of distance from 
the generating point did not permit identification of the interface 
phases . 

Neither model is by any means limited to the homogeneous layer- 
transition zone structure presented in this paper. The velocity 
structure could consist entirely of transition regions of various 
gradients. This could be applicable to investigation of the effects of 
the gradients on reverberation. 

A collection of records is on file at the Naval Postgraduate School, 
Department of Oceanography (Code 58Ad) . 
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KK SUCCESSIVE GRAPHS OF STRESS VERSUS TIME 
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NUMBER OF POINTS DEFINING THE WAVE VELOCITY AND DENSITY 
POLYGONS, NC,LE,?0 
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SUBROUTINE GRin( XMINt XMAX, YM I N , YMA X , GX ,G Y , X , Y , KA , L ) 

THIS SUBROUTINE GENERATES A SEQUENCE OF POINTS (X,Y) 
WHICH MAY BE USEO TO PLOT A GRID ON A CALCQMP PLOTTER 
GRAPH 
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NOTE: A 6 BY 10 INCH GRID OF ONE INCH SQUARES REQUIRES 

36 POINTS. 

AN 8 BY 60 INCH GRID OF HALF INCH SQUARES REQUIRES 
280 POINTS. 
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